#!/usr/bin/awk -f
#################################################################
#
# makeapot:  generate analytic potential files from scratch
#
#################################################################
#
#   Copyright 2009-2010 Daniel Schopf
#             Institute for Theoretical and Applied Physics
#             University of Stuttgart, D-70550 Stuttgart, Germany
#             http://www.itap.physik.uni-stuttgart.de/
#
#################################################################
#
#   This file is part of potfit.
#
#   potfit is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   potfit is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with potfit; if not, see <http://www.gnu.org/licenses/>.
#
#################################################################
#
# Usage: makeapot <options>
# 	 -cp # 			enable chemical potential
# 	 -co # 			cutoff
# 	 -n # type 		number of potentials of this type
# 				(_sc can be used)
#
# eg:    makeapot -cp -cn 5 MgZn.config -co 5 -n 3 eopp_sc
#
#################################################################

function write_header(n)
{
    if (n==1 || n==6)
	    inter="PAIR";
    else if (n==7 || n==12)
	    inter="EAM";
    else if (n==5 || n==13 || n==24)
	    inter="ADP";
    else if (n==3) {
	print "potential type is not unique for 3 potentials, assuming EAM" > "/dev/stderr";
	inter="EAM";
	    }
    else {
	printf "Unknown number (%d) of potential functions, aborting.\n",n > "/dev/stderr";
        exit;
	    }
    print "#F 0 "n;
    printf "#T %s\n",inter;
    printf "#I";
    for (i=0;i<n;i++)
	printf " 0";
    printf "\n#E\n\n";
}

function write_cp(n)
{
    for (i=0;i<n;i++)
	printf "cp_%d -1 -10 0\n",i;
    printf "\n";
}

function write_cn(n)
{
    print "cn "n;
    if (!system("[ -e "cn_file" ]")) {
	("compnodes -n "n" "cn_file) | getline x;
	split(x,mu);
    } else {
	for (i=0;i<n;i++)
	    mu[i+1]=(1/n)*(i+1);
    }
    for (i=0;i<n;i++)
	printf "%.3f -1 -10 0\n",mu[i+1];
}

function valminmax(val,min,max)
{
  if (random)
	temp = sprintf("%.2f %.2f %.2f\n",min+(max-min)*rand(),min,max);
  else
	temp = sprintf("%.2f %.2f %.2f\n",val,min,max);
  return temp;
}

function write_pot(name)
{
    print "type "name;
    print "cutoff "co;
    do_smooth=0;
    if (match(name,"_sc$")>0) {
	gsub("_sc$","",name);
	do_smooth=1;
    }
    if (name=="eopp") {
	printf "C_1 %s",valminmax(15,1,10000);
	printf "eta_1 %s",valminmax(6,1,20);
	printf "C_2 %s",valminmax(5,-100,100);
	printf "eta_2 %s",valminmax(3,1,10);
	printf "k %s",valminmax(2.5,0,6);
	printf "phi %s",valminmax(3,0,6.3);
    }
    else if (name=="morse") {
	printf "D %s",valminmax(0.1,0,1);
	printf "alpha %s",valminmax(2,1,5);
	printf "r_0 %s",valminmax(2.5,1,5);
    }
    else if (name=="power_decay") {
	printf "a %s",valminmax(1,0.1,10);
	printf "b %s",valminmax(2,1,5);
    }
    else if (name=="bjs") {
	printf "F_0 %s",valminmax(-1,-10,0);
	printf "gamma %s",valminmax(2,0.1,2);
	printf "F_1 %s",valminmax(0,-10,10);
    }
    else if (name=="parabola") {
	printf "alpha %s",valminmax(1,-10,10);
	printf "beta %s",valminmax(1,-10,10);
	printf "gamma %s",valminmax(1,-10,10);
    }
    else if (name=="csw") {
	printf "a_1 %s",valminmax(0.2,-2,2);
	printf "a_2 %s",valminmax(0.2,-2,2);
	printf "alpha %s",valminmax(2,1,6);
	printf "beta %s",valminmax(3,0.5,5);
    }
    else if (name=="csw2") {
	printf "a %s",valminmax(0.2,-2,2);
	printf "alpha %s",valminmax(2,1,6);
	printf "phi %s",valminmax(0,0,6.3);
	printf "beta %s",valminmax(3,0.5,5);
    }
    else if (name=="universal") {
	printf "F_0 %s",valminmax(-1,-10,10);
	printf "m %s",valminmax(1,0,20);
	printf "n %s",valminmax(2,0,20);
	printf "F_1 %s",valminmax(0,0,0);
    }
    else if (name=="lj") {
	printf "epsilon %s",valminmax(0.1,0,1);
	printf "sigma %s",valminmax(2.5,1,4);
    }
    else if (name=="softshell") {
	printf "alpha %s",valminmax(1,0.1,10);
	printf "beta %s",valminmax(2,1,5);
    }
    else if (name=="exp_decay") {
	printf "alpha %s",valminmax(4,1,20);
	printf "beta %s",valminmax(1,0.5,5);
    }
    else if (name=="eopp_exp") {
	printf "C_1 %s",valminmax(15,0.5,10000);
	printf "eta_1 %s",valminmax(6,1,20);
	printf "C_2 %s",valminmax(5,-100,100);
	printf "eta_2 %s",valminmax(3,1,20);
	printf "k %s",valminmax(2.5,0,6);
	printf "phi %s",valminmax(3,0,6.3);
    }
    else if (name=="meopp") {
	printf "C_1 %s",valminmax(15,0.5,10000);
	printf "eta_1 %s",valminmax(6,1,20);
	printf "C_2 %s",valminmax(5,-100,100);
	printf "eta_2 %s",valminmax(3,1,20);
	printf "k %s",valminmax(2.5,0,6);
	printf "phi %s",valminmax(3,0,6.3);
	printf "r_0 %s",valminmax(0,-3,3);
    }
    else if (name=="const") {
	printf "c %s",valminmax(1,0,10);
    }
    else if (name=="poly_5") {
	printf "F_0 %s",valminmax(1,-100,100);
	printf "F_1 %s",valminmax(1,-100,100);
	printf "q_1 %s",valminmax(1,-100,100);
	printf "q_2 %s",valminmax(1,-100,100);
	printf "q_3 %s",valminmax(1,-100,100);
    }
    else if (name=="double_morse") {
	printf "E_1 %s",valminmax(1,-10,10);
	printf "alpha_1 %s",valminmax(1,-10,10);
	printf "r_0 %s",valminmax(1,-10,10);
	printf "E_2 %s",valminmax(1,-10,10);
	printf "alpha_2 %s",valminmax(1,-10,10);
	printf "r_1 %s",valminmax(1,-10,10);
	printf "delta %s",valminmax(1,-10,10);
    }
    else if (name=="double_exp") {
	printf "a %s",valminmax(1,-10,10);
	printf "beta_1 %s",valminmax(1,-10,10);
	printf "r_0 %s",valminmax(1,-10,10);
	printf "beta_2 %s",valminmax(1,-10,10);
	printf "r_1 %s",valminmax(1,-10,10);
    }
    else if (name=="cbb") {
        printf "z %s",valminmax(3,-6,6);
	printf "f0b %s",valminmax(1.5,0,3);
	printf "a %s",valminmax(3,2,4);
	printf "b %s",valminmax(0.5,0,1);
	printf "c %s",valminmax(300,200,400);
	printf "f0d %s",valminmax(50,40,60);
	printf "beta %s",valminmax(1.6,1.4,1.8);
	printf "r_0 %s",valminmax(2.5,2.25,2.75);
    }
    else if (name=="exp_plus") {
        printf "d1 %s",valminmax(1,0,10);
	printf "d2 %s",valminmax(1,-10,10);
	printf "d3 %s",valminmax(0,-2,2);
    }
    if (do_smooth)
        if (global_co)
	    printf "h!\n";
        else
	    printf "h %s",valminmax(1,0.5,2);
}

BEGIN{
    co = 7;
    j = 0;
    total_pot = 0;
    for (i=1;i<ARGC;i++) {
	if (ARGV[i]=="-cp")
	    enable_cp=1;
	if (ARGV[i]=="-h")
	    global_co=1;
	if (ARGV[i]=="-cn") {
	    enable_cn=1;
	    cn=ARGV[i+1];
	    cn_file=ARGV[i+2];
	}
	if (ARGV[i]=="-n") {
	    total_pot+=ARGV[i+1];
	    for (k=0;k<ARGV[i+1];k++){
		pot[j]=ARGV[i+2];
		j++;
	    }
	}
	if (ARGV[i]=="-co")
	    co=ARGV[i+1];
	if (ARGV[i]=="-r")
	    random=1;
    }
    srand(strftime( "%s", systime() ));
    write_header(total_pot);
    if (enable_cp) {
	write_cp(-0.5+sqrt(0.25+2*total_pot));
	if (enable_cn)
	    write_cn(cn);
    }
    if (global_co) {
	printf "global 1\n";
	printf "h %s",valminmax(1,0.5,2);
	printf "\n";
	    }
    for (i=0;i<total_pot;i++) {
	if (pot[i]=="pohlong")
		pot[i]="bjs";
	write_pot(pot[i]);
	if (i!=(total_pot-1))
		printf "\n";
    }

}

